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Abstract. This review article provides an overview of recent work in 
the modeling and analysis of recurrent events arising in engineering, 
reliability, public health, biomedicine and other areas. Recurrent event 
modeling possesses unique facets making it different and more diffi- 
cult to handle than single event settings. For instance, the impact of 
an increasing number of event occurrences needs to be taken into ac- 
count, the effects of covariates should be considered, potential asso- 
ciation among the interevent times within a unit cannot be ignored, 
and the effects of performed interventions after each event occurrence 
need to be factored in. A recent general class of models for recurrent 
events which simultaneously accommodates these aspects is described. 
Statistical inference methods for this class of models are presented and 
illustrated through applications to real data sets. Some existing open 
research problems are described. 

Key words and phrases: Counting process, hazard function, martin- 
gale, partial likelihood, generalized Nelson-Aalen estimator, general- 
ized product-limit estimator. 



1. INTRODUCTION 

A decade ago in a Statistical Science article, 
Singpurwalla (1995) advocated the development, 
adoption and exploration of dynamic models in the 
theory and practice of reliability. He also pinpointed 
that the use of stochastic processes in the model- 
ing of component and system failure times offers a 
rich environment to meaningfully capture dynamic 
operating conditions. In this article, we review re- 
cent research developments in dynamic failure-time 
models, in the context of both stochastic model- 
ing and inference concerning model parameters. Dy- 
namic models are not limited in applicability and 
relevance to the engineering and reliability areas. 
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They are also relevant in other fields such as public 
health, biomedicine, economics, sociology and poli- 
tics. This is because in many studies in these var- 
ied areas, it is oftentimes of interest to monitor the 
occurrences of an event. Such an event could be the 
malfunctioning of a mechanical or electronic system, 
encountering a bug in computer software, the out- 
break of a disease, occurrence of a migraine, reoc- 
currence of a tumor, a drop of 200 points in the Dow 
Jones Industrial Average during a trading day, com- 
mission of a criminal act in a city, serious disagree- 
ments between a married couple, or the replacement 
of a cabinet minister /secretary in an administration. 
These events recur and so it is of interest to de- 
scribe their recurrence behavior through a stochastic 
model. If such a model has excellent predictive abil- 
ity for event occurrences, then if event occurrences 
lead to drastic and/or negative consequences, pre- 
ventive interventions may be attempted to minimize 
damages, whereas if they lead to beneficial and/or 
positive outcomes, certain actions may be performed 
to hasten event occurrences. 

It is because of performed interventions after event 
occurrences that dynamic models become especially 
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pertinent, since through such interventions, or some- 
times noninterventions, the stochastic structure gov- 
erning future event occurrences is altered. This change 
in the mechanism governing the event occurrences 
could also be due to the induced change in the struc- 
ture function in a reliability setting governing the 
system of interest arising from an event occurrence 
(cf. Hollander and Peha, 1995; Kvam and Peha, 2005; 
Peha and Slate, 2005). Furthermore, since several 
events may occur within a unit, it is also important 
to consider the association among the interevent 
times which may arise because of unobserved ran- 
dom effects or frailties for the unit. In addition, 
there is a need to take into account the potential 
impact of environmental and other relevant covari- 
ates, possibly including the accumulated number of 
event occurrences, which could affect future event 
occurrences. 

In this review article we describe a flexible and 
general class of dynamic stochastic models for event 
occurrences. This class of models was proposed by 
Peha and Hollander (2004). We also discuss infer- 
ence methods for this class of models as developed 
in Peha, Strawderman and Hollander (2001), Peha, 
Slate and Gonzalez (2007) and Kvam and Peha (2005). 
We demonstrate their applications to real data sets 
and indicate some open research problems. 

2. BACKGROUND AND NOTATION 

Rapid and general progress in event time mod- 
eling and inference has benefited immensely from 
the adoption of the framework of counting processes, 
martingales and stochastic integration as introduced 
in Aalen's (1978) pioneering work. The present re- 
view article adopts this mathematical framework. 
Excellent references for this framework are the mono- 
graph of Gill (1980), and the books by Fleming and 
Harrington (1991) and Andersen, Borgan, Gill and 
Keiding (1993). We introduce in this section some 
notation and very minimal background in order to 
help the reader in the sequel, but advise the reader 
to consult the aforementioned references to gain an 
in-depth knowledge of this framework. 

We denote by (0, J 7 , P) the basic probability space 
on which all random entities are defined. Since in- 
terest will be on times between event occurrences 
or on the times of event occurrences, nonnegative- 
valued random variables will be of concern. For a 
random variable T with range = [0, oo), F(t) = 
F T (t) = P{T < t} and S(t) = S T (t) = 1 - F(t) = 



P{T > t} will denote its distribution and survivor 
(also called reliability) functions, respectively. The 
indicator function of the event A will be denoted 
by /{A}. The cumulative hazard function of T is 
defined according to 



A(t)=A T (t) = I{t>0} f 

Jo 



F(dw) 



S(w-) 

with the convention that S(w—) = lim t + w S(t) = 
P{T > w} and /„••• = J( 6] " " ■ ^ or a non decreasing 
function G:$l + — ► 3? + with G(0) = 0, its product- 
integral over (0, t] is defined via (cf. Gill and Jo- 
hansen, 1990; Andersen et al., 1993) H^ll-Gidw)) 

hm M ^ooII=i[l " (G(U) ~ G(ti_i))], where as 
M — > oo, the partition = to < t\ < ■ ■ ■ < tj[f = to 
satisfies maxi<j<jvf \U — tj— 1| — > 0. The survivor func- 
tion in terms of the cumulative hazard function be- 
comes 

t 

(2.1) S(t) = I{t < 0} + I{t > 0} {J [1 - A(diu)]. 

w=0 

For a discrete random variable T with jump points 
{tj}, the hazard Xj at tj is the conditional probabil- 
ity that T = tj, given T >tj, so A(t) = J2{j ■ tj <t} -\? • 
If T is an absolutely continuous random variable 
with density function f(t), its hazard rate function is 
A(t) = f(t)/S(t), so A(t) =J* \(w)dw = -logS(t). 
The product-integral representation of S(t) accord- 
ing to whether T is purely discrete or purely contin- 
uous is 



S(t)=U[l-A(dw)} 
(2.2) 

IB 1 -Ail. 

tn<t 



if T is discrete, 
exp| — J \(w)dvu^, if T is continuous. 



A benefit of using hazards or hazard rate functions 
in modeling is that they provide qualitative aspects 
of the event occurrence process as time progresses. 
For instance, if the hazard rate function is increas- 
ing (decreasing) then this indicates that the event 
occurrences are increasing (decreasing) as time in- 
creases, and thus we have the notion of increasing 
(decreasing) failure rate [IFR (DFR)] distributions. 
For many years, it was the focus of theoretical re- 
liability research to deal with classes of distribu- 
tions such as IFR, DFR, increasing (decreasing) fail- 
ure rate on average [IFRA (DFRA)], and so on, 
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specifically with regard to their closure properties 
under certain reliability operations (cf. Barlow and 
Proschan, 1981). 

In monitoring an experimental unit for the oc- 
currence of a recurrent event, it is convenient and 
advantageous to utilize a counting process {N(s), 
s > 0}, where N(s) denotes the number of times that 
the event has occurred over the interval [0, s] . The 
paths of this stochastic process are step-functions 
with iV(0) = and with jumps of unity. If we rep- 
resent the calendar times of event occurrences by 
Si < S2 < £3 < ■ ■ ■ with the convention that So = 0, 
then N(s) = Y^,j? =1 I{Sj The interevent times 
are denoted by Tj = Sj — Sj-i, j = 1, 2, In speci- 
fying the stochastic characteristics of the event 
occurrence process, one either specifies all the finite- 
dimensional distributions of the process {iV(s)}, or 
specifies the joint distributions of the Sj's or the 
Tj's. For example, a common specification for event 
occurrences is the assumption of a homogeneous Pois- 
son process (HPP) where the interevent times Tj are 
independent and identically distributed exponential 
random variables with scale (3. This is equivalent 
to specifying that, for any s\ < S2 < ■ • • < sk, the 
random vectors (N(si), N(s2) — N(si), . . . , N(sk) — 
N(sk-i)) have independent components with N(sj) — 
N(sj—i) having a Poisson distribution with param- 
eter (3{sj — Sj-i). From this specification, the finite- 
dimensional distributions of (N(si), AT(s 2 ), • • • , 
N(sk)) can be obtained. 

An important concept in dynamic modeling is that 
of a history or a filtration, a family of <r-fields F = 
{J~s '■ s > 0} where J- s is a sub-u-field of J- with T Sl C 
T S2 for every s\ < s 2 and with T s = f] h ^ F s +h, a 
right-continuity property. One interprets J- s as all 
information that accrued over [0, s] . When augmented 
in (O, P), we obtain the filtered probability space 
(£l,J-,F,P). It is with respect to such a filtered 
probability space that a process {X(s) : s > 0} is said 
to be adapted [X(s) is measurable with respect to 
•7"sj Vs > 0]; is said to be a (sub)martingale [adapted, 
E\X(s)\ < 00, and, Vs,i > 0, E{X(s + t)\T s }(>) = 
X(s) almost surely]. Doob-Meyer's decomposition 
guarantees that for a submartingale Y = {Y(s) : s > 
0} there is a unique increasing predictable (measur- 
able with respect to the u-field generated by adapted 
processes with left-continuous paths) process A = 
{A(s) : s > 0}, called the compensator, with A(0) = 
such that M = {M(s) = Y(s) - A(s):s > 0} is 
a martingale. Via Jensen's inequality, then for a 
square-integrable martingale X = {X(s):s > 0}, 



there is a unique compensator process A = {^4(s) : s > 
0} such that X 2 — A is a martingale. This process A, 
denoted by (X) = {{X)(s) : s > 0}, is called the pre- 
dictable quadratic variation (PQV) process of X. A 
useful heuristic, though somewhat imprecise, way of 
presenting the main properties of martingales and 
PQV's is through the following differential forms. 
For a martingale M, E{dM{s)\T s -} = 0,Vs > 0; 
whereas for the PQV (M), E{dM 2 {s)\T s ^} = 
\&i{dM{s)\r s ^} = d{M)(s)ys > 0. 

For the HPP N = {N(s) : s > 0} with rate (3 and 
with F = {Ts = a{N{w),w < s} : s > 0}, N is a sub- 
martingale since its paths are nondecreasing. Its com- 
pensator process is A = {A(s) = (3s : s > 0}, so that 
M = {M(s) = N(s) -/3s:s>0} is a martingale. 
Furthermore, since N(s) is Poisson-distributed with 
rate (3s, so that {M 2 (s) - A{s) = (N(s) - (3s) 2 - 
(3s : s > 0} is a martingale, the PQV process of M 
is also A. Through the heuristic forms, we have 
E{dN(s)\F s -} = dA(s),s > 0. Since dN(s) G {0, 1}, 
then we obtain the probabilistic expression P{dN(s) = 
IIj^} = dA(s),s > 0. Analogously, Var{diV(s)j 
Ts-} = E{[dN(s) - dAis)} 2 ^} = dA(s),s > 0. 
These formulas hold for a general counting process 
{N(s) : s > 0} with compensator process {A(s) : s > 
0}. In essence, conditionally on J- s -, dN(s) has a 
Bernoulli distribution with success probability dA(s). 
Over the interval [0, s] , following Jacod, the like- 
lihood function can be written in product-integral 
form as 

s 

L(s) = [] {dA(w)} dN{w) {l - dA(w)} 1 - dIf W 
= \f[[dA(w)] dN ^\exp{-A(s)}, 

lw=0 ) 

with the last equality holding when A(-) has contin- 
uous paths. 

Stochastic integrals play a crucial role in this stochas- 
tic process framework for event time modeling. For 
a square-integrable martingale X = {X(s):s > 0} 
with PQV process (X) = {(X)(s):s> 0}, and for 
a bounded predictable process Y = {Y(s):s > 0}, 
the stochastic integral of Y with respect to X , de- 
noted by / YdX = {/o Y(w) dX{w) : s > 0}, is well 
defined. It is also a square-integrable martingale with 
PQV process (J YdX) = {Jq Y 2 (w) d(X)(w) : s > 0}. 
When X is associated with a counting process N, 
that is, X = N — A, the paths of the stochastic in- 
tegral jYdX can be taken as pathwise Lebesgue- 
Stieltjes integrals. 
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Martingale theory also plays a major role in ob- 
taining asymptotic properties of estimators as first 
demonstrated in Aalen (1978), Gill (1980) and An- 
dersen and Gill (1982). The main tools used in asymp- 
totic analysis are Lenglart's inequality (cf. Lenglart, 
1977; Fleming and Harrington, 1991; Andersen et 
al., 1993) which is used in proving consistency, and 
Rebolledo's (1980) martingale central limit theorem 
(MCLT) (cf. Fleming and Harrington, 1991; Ander- 
sen et al., 1993) which is used for obtaining weak 
convergence results. We refer the reader to Fleming 
and Harrington (1991) and Andersen et al. (1993) 
for the in-depth theory and applications of these 
modern tools in failure-time analysis. 

3. CLASS OF DYNAMIC MODELS 

Let us now consider a unit being monitored over 
time for the occurrence of a recurrent event. The 
monitoring period could be a fixed interval or it 
could be a random interval, and for our purposes 
we denote this period by [0,r], where r is assumed 
to have some distribution G, which may be degen- 
erate. With a slight notational change from Sec- 
tion 2 we denote by N^(s) the number of events 
that have occurred on or before time s, and by Y' (s) 
an indicator process which equals 1 when the unit 
is still under observation at time s, otherwise. 
With S*o = and Sj, j = 1,2, ... , denoting the suc- 
cessive calendar times of event occurrences, and with 
Tj = Sj — Sj-\,j = 1,2,..., being the interevent or 
gap times, observe that 

oo 

N\s) = J2l{Sj <min(s,r)} and 



Associated with the unit is a, possibly time-depen- 
dent, lxg covariate vector X = {X(s) : s > 0}. In re- 
liability engineering studies, the components of this 
covariate vector could be related to environmental 
or operating condition characteristics; in biomedical 
studies, they could be blood pressure, treatment as- 
signed, initial tumor size, and so on; in a sociological 
study of marital disharmony, they could be length 
of marriage, family income level, number of children 
residing with the couple, ages of children, and so on. 
Usually, after each event occurrence, some form of 
intervention is applied or performed, such as replac- 
ing or repairing failed components in a reliability 
system, or reducing or increasing physical activity 



after a heart attack in a medical setting. These inter- 
ventions will typically impact the next occurrence of 
the event. There is furthermore recognition that for 
the unit the interevent times are associated or cor- 
related, possibly because of unobserved random ef- 
fects or so-called frailties. A pictorial representation 
of these aspects is given in Figure 1. Observe that 
because of the finiteness of the monitoring period, 
which leads to a sum-quota accrual scheme, there 
will always be a right-censored interevent time. The 
observed number of event occurrences over [0,r], 
K = N\t), is also informative about the stochas- 
tic mechanism governing event occurrences. In fact, 
since K = maxjfc : Y^j=i Tj < t}, then, conditionally 
on (K, r), the vector (Ti, Ti, ■ ■ ■ , Tk) has dependent 
components, even if at the outset TijTJj,... are in- 
dependent. 

Recognizing these different aspects in recurrent 
event settings, Peha and Hollander (2004) proposed 
a general class of models that simultaneously incor- 
porates all of these aspects. To describe this class 
of models, we suppose that there is a filtration F = 
{Fs-s > 0} such that for each s > 0, a{(N^(v), 
Y^(v+),X(v+),£(v+)):v <s}C JT S . We also assume 
that there exists an unobservable positive random 
variable Z, called a frailty, which induces the cor- 
relation among the inter-event times. The class of 
models of Peha and Hollander (2004) can now be 
described in differential form via 

P{dN(s) = l\F s -,Z} 

(3.5) = ZY^s)X Q [£(s)} 

■ p[Ni(s-);a]rp(X(s)0)ds, 

where Ao(-) is a baseline hazard rate function; p(-; •) 
is a nonnegative function with p(0; •) = 1 and with 
a being some parameter; ip(-) is a nonnegative link 
function with f5 a q x 1 regression parameter vector; 
and Z is a frailty variable. The at-risk process Y*(s) 
indicates that the conditional probability of an event 
occurring becomes zero whenever the unit is not un- 
der observation. Possible choices of the />(■;•) and 
Tp(-) functions are p(k;a) = a k and ip{w) =exp(u>), 
respectively. For the geometric choice of p(-;-), if 
a > 1 the effect of accumulating event occurrences 
is to accelerate event occurrences, whereas if a < 1 
the event occurrences decelerate, the latter situa- 
tion appropriate in software debugging. The process 
£{■) appearing as argument in the baseline hazard 
rate function, called the effective age process, is pre- 
dictable, observable, nonnegative and pathwise al- 
most surely differentiable with derivative £'(■)■ This 
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effective age process models the impact of performed 
interventions after each event occurrence. A picto- 
rial depiction of this effective age process is in Fig- 
ure 2. In this plot, after the first event, the per- 
formed intervention has the effect of reverting the 
unit to an effective age equal to the age just be- 



fore the event occurrence (called a minimal repair 
or intervention), while after the second event, the 
performed intervention has the effect of reverting 
the effective age to that of a new unit (hence this 
is called a perfect intervention or repair). After the 
third event, the intervention is neither minimal nor 




Covariate vector: X(s) - (X^s), . . . , Xg(s)) 

Fig. 1. Pictorial depiction of the recurrent event data accrual for a unit illustrating the window of observation [0, t], inter- 
vention performed after an event occurrence, an unobserved frailty Z , the presence of a vector of covariates X, the interevent 
times Tj and the calendar times of event occurrences Sj . 




Calendar Time L 

Fig. 2. An example of an effective age process, S(s), for a unit encountering successive occurrences of a recurrent event. 
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perfect and it has the effect of restarting the effective 
age at a value between zero and that just before the 
event occurred; while for the fourth event, the inter- 
vention is detrimental in that the restarting effective 
age exceeds that just before the event occurred. 

The effective age process could occur in many 
forms, and the idea is this should be determined 
in a dynamic fashion in conjunction with interven- 
tions that are performed. As such the determination 
of this process should preferably be through con- 
sultations with experts of the subject matter under 
consideration. Common forms of this effective age 
process are: 

Minimal Intervention or Repair: 

(3.6) 

£(s) = s, 
Perfect Intervention or Repair: 

(3.7) 

£(*) = s-5jvt( s -), 
BBS Model: 

(3.8) 

£(s) = s-Sr v{s _ ) , 

where in (3.8) with /i,l2,... being independent 
Bernoulli random variables with 1^ having success 
probability p(Sk) withp: 9? + — ► [0, 1], we define r/(s) = 

Eili^ h an d r = 0, T k = min-fj > : 

Ij = l},k = 1,2, Thus, in (3.8), £(s) represents 

the time measured from s since the last perfect re- 
pair. This effective age is from Block, Borges and 
Savits (1985), whereas when p(s) =p for some p G 
[0, 1] , we obtain the effective age process for the 
Brown and Proschan (1983) minimal repair model. 
Clearly, the effective age functions (3.6) and (3.7) 
are special cases of (3.8). Other effective age pro- 
cesses that could be utilized are those associated 
with the general repair model of Kijima (1989), Stadje 
and Zuckerman (1991), Baxter, Kijima and Tortorella 
(1996), Dorado, Hollander and Sethuraman (1997) 
and Last and Szekli (1998). See also Lindqvist (1999) 
and Lindqvist, Elvebakk and Heggland (2003) for 
a review of some of these models pertaining to re- 
pairable systems, and Gonzalez, Peha and Slate (2005) 
for an effective age process suitable for cancer stud- 
ies. 

A crucial property arising from the intensity spec- 
ification in (3.5) is it amounts to postulating that, 
with 

A\s;Z,\ (-),a,(3) 



(3.9) = Z f Y\w)\ Q [£(w)\ 
Jo 

• p[N*(w-);a]il>(X.(w)P) dw, 
then, conditionally on Z, the process 
{Mt( S ;Z,A o (0,a,/3) 

(3.10) 

= N^s)-A\s;Z,X (-),a,(3):s>0} 

is a square-integrable martingale with PQV process 
A'(-; Z, Ao(-), a, [3). The class of models is general 
and flexible and subsumes many current models for 
recurrent events utilized in survival analysis and reli- 
ability. In particular, it includes those of Jelinski and 
Moranda (1972), Gail, Santner and Brown (1980), 
Gill (1981), Prentice, Williams and Peterson (1981), 
Lawless (1987), Aalen and Husebye (1991), Wang 
and Chang (1999), Peha, Strawderman and Hollan- 
der (2001) and Kvam and Peha (2005). We demon- 
strate the class of models via some concrete exam- 
ples. 

Example 1 . The first example concerns a load- 
sharing i^-component parallel system with identi- 
cal components. The recurrent event of interest is 
component failure and failed components are not re- 
placed. When a component fails, a redistribution of 
the system load occurs among the remaining func- 
tioning components, and to model this system in 
a general way, we let otQ = I and a±, . . . ,ok-i be 
positive constants, referred to as load-share param- 
eters. One possible specification of these parameters 
is a k = K/ (K - k) , k = 0, 1 , 2, . . . , K - 1 , though they 
could be unknown constants, possibly ordered. The 
hazard rate of event occurrence at calendar time s, 
provided that the system is still under observation, 
is A(s) = Xq(s)[K - A^ t (s-)]a Art ( s -), where A (-) is 
the common hazard rate function of each compo- 
nent and A^(s) is the number of component fail- 
ures observed on or before time s. This is a special 
case of the general class of models with £ (s) = s, 
p(k; a±, . . . , OiK-i) = [K — k]a k , and ijj(w) = 1. This 
is the equal load-sharing model in Kvam and Peha 
(2005). More generally, this could accommodate en- 
vironmental or operating condition covariates for 
the system, and even an unobserved frailty compo- 
nent. 

Example 2. Assume in a software reliability model 
that there are a bugs at the beginning of a debug- 
ging process and the event of interest is encounter- 
ing a bug. A possible model is these a bugs are 
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competing to be encountered, and if each of them 
has hazard rate of Ao(s) of being encountered at 
time s, then the total hazard rate at time s of the 
software failing is Xo(s)a. Upon encountering a bug 
at time Si, this bug is eliminated, thus decreasing 
the number of remaining bugs by one. The debug- 
ging process is then restarted at time just after Si 
(assuming it takes zero time to eliminate the bug, 
clearly an oversimplification). In general, suppose 
that just before calendar time s, N^(s-) bugs have 
been removed, and the last bug was encountered 



at calendar time S 



ATt(s-) 



. Then, the overall hazard 



of encountering a bug at calendar time s with s > 
^jvt( s -) is \ (s - S N t( s -))[a - N\s-)\. Thus, pro- 
vided that the debugging process is still in progress 
at time s, then the hazard of encountering a bug at 
time s is A(s) = Xq{s — S , Ar t( s _)) max[0, a — iV' (s— )]. 
Again, this is a special case of the general class 
of models with £(s) = s — S N u s -), a consequence 
of the restart of the debugging process, p(k; a) = 
max{0,a — k} and ip(w) = 1. This software debug- 
ging model is the model of Jelinski and Moranda 
(1972) and it was also proposed by Gail, Santner 
and Brown (1980) as a carcinogenesis model. See 
also Agustin and Peha (1999) for another model in 
software debugging which is a special case of the 
general class of models. 

Cox's (1972) proportional hazards model is one 
of the most used models in biomedical and public 
health settings. Extensions of this model have been 
used in recurrent event settings, and Therneau and 
Grambsch (2000) discuss some of these Cox-based 
models such as the independent increments model 
of Andersen and Gill (1982), the conditional model 
of Prentice, Williams and Peterson (1981) and the 
marginal model of Wei, Lin and Weissfeld (1989). 
The independent increments model is a special case 
of the general class of models obtained by taking ei- 
ther £(s) =s or £(s) = s — S^u s _\ with p{k;a) = 
1 and ip(w) = exp(u;). The marginal model strati- 
fies according to the event number and assumes a 
Cox-type model for each of these strata, with the 
jth interevent time in the ith unit having intensity 
Yij(t)Xoj(t) exp{Xi(t)Pj} , where Yij(t) equals 1 un- 
til the occurrence of the jth event or when the unit 
is censored. The conditional model is similar to the 
marginal model except that Yij{t) becomes 1 only 
after the (j — l)st event has occurred. 



4. STATISTICAL INFERENCE 

The relevant parameters for the model in (3.5) are 
Ao(-)> a, /3 and the parameter associated with the 
distribution of the frailty variable Z. A variety of 
forms for this frailty distribution is possible, but we 
restrict to the case where Z has a gamma distribu- 
tion with mean 1 and variance l/£. The parame- 
ter associated with G, the distribution of r, is usu- 
ally viewed as a nuisance, though in current joint 
research with Akim Adekpedjou, a Ph.D. student 
at the University of South Carolina, the situation 
where G is informative about the distributions of 
the interevent times is being explored. 

Knowing the values of the model parameters is im- 
portant because the model can be utilized to predict 
future occurrences of the event, an important issue 
especially if an event occurrence is detrimental. To 
gain knowledge about these parameters, a study is 
performed to produce sample data which is the ba- 
sis of inference about the parameters. We consider a 
study where n independent units are observed and 
with the observables over (calendar) time [0, s*] de- 
noted by 

DATA n (s*) 

(4.11) ={[(X i (v),Nj(v),Y^v),S i (v)),v<s*], 

i = l,2,..., n}. 

Observe that DATA n (s*) provides information 
about the Tj's. More generally, we observe the ni- 
trations {(Fiv, v < s*), i = 1, 2, . . . , n} or the overall 
filtration F s » = {T v = \f? = iF iv ,v <s*}. The goals 
of statistical inference are to obtain point or inter- 
val estimates and/or test hypotheses about model 
parameters, as well as to predict the time of occur- 
rence of a future event, when DATA n (s*) or F s * is 
available. We focus on the estimation problem be- 
low, though we note that the prediction problem is 
of great practical importance. 

Conditional on Z = [Zi, Z2, ■ ■ ■ , Z n ), from (2.3) the 
likelihood process for (Ao(-)>«> 0) is 



(4.12) 



L c {s;Z,X (-),a,f3) 

■-f[{z^lf[m V] x (-), a , m dN ^ 



i=l 



.v=0 



■ exp 



Z, 



n 



Bi(v;X (-),a,(3) dv 



where B^v; X (-),a, f3) = Yf(v)\o[£i(v)]p[N}(v-); 
a]?/4Xj(f )0\. Observe that the likelihood process when 
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the model does not involve any frailties is obtained 
from (4.12) by setting = 1, i = 1, 2, . . . , n, which is 
equivalent to letting £ — > oo. The resulting no- frailty 
likelihood process is 

L(s;\ (-),a,(3) 
( 4 - 13 ) =f[{{f[mv;Xo(-U,p)] dNl M 



.t)=0 



• exp 



Bi(v;\o(-),a,{3) dv 



This likelihood process is the basis of inference about 
(Xo(-),oc,P) in the absence of frailties. Going back 
to (4.12), by marginalizing over Z under the gamma 
frailty assumption, the likelihood process for (Ao(-)> a, 
becomes 

L(s;X (-),a,(3,£) 




i=l 



£ + Jo ^(^S A o(-)> a > Z 3 ) rfw 



n 

D = 



l-idJvf(v)' 



(4.14) 



£ + y^ Bi(w;X (-),a,j3)dw 

There are two possible specifications for the base- 
line hazard rate function Ao(-) : parametric or non- 
parametric. If parametrically specified, then it is 
postulated to belong to some parametric family of 
hazard rate functions, such as the Weibull or gamma 
family, that depends on an unknown p x 1 param- 
eter vector 0. In this situation, a possible estima- 
tor of (0, a, /?,£), given F s *, is the maximum like- 
lihood estimator (0, a,/3, £), the maximizer of the 
mapping (0, a, (3, £) i— » L(s*; Ao(-; 0), oc,(3, £)• Stocker 
(2004) studied the finite- and large-sample prop- 
erties, and associated computational issues, of this 
parametric ML estimator in his dissertation research. 
In particular, following the approach of Nielsen, Gill, 
Andersen and S0rensen (1992), he implemented an 
expectation-maximization (EM) algorithm (cf. 
Dempster, Laird and Rubin, 1977) for obtaining the 
ML estimate. In this EM implementation the frailty 
variates Z{ are viewed as missing and a variant of 
the no- frailty likelihood process in (4.13) is used 
for the maximization step in this algorithm. We re- 
fer to Stocker (2004) for details of this computa- 
tional implementation. For large n, and under cer- 
tain regularity conditions, it can also be shown that 



(0, is approximately multivariate normally 
distributed with mean (0, a, (3, £) and covariance ma- 
trix i/ _1 (#,a,/3,£), where I(0,a,f3,£) is the ob- 
served Fisher information associated with the like- 
lihood function L(s*; Ao(-; 0), a, /?, £). That is, with 

e = (e, a ,M\ i(0, a, (3,o = -{p/deae*} ■ 

l(s*;Xo(-;0),a,(3,£), where l(s; Ao(-; 0),a,(3,£,) is the 
log-likelihood process given by 



l(s;X o (-,0),a,(3,O 



(4.15) + 



i + Jo B i( w \ A o(-5 0), a, (3) dw 



log 



(N}(v-) + 



£+ / Bi(w } X o (-,0),a,f3)dw 



dN}(i 



Tests of hypotheses and construction of confidence 
intervals about model parameters can be developed 
using the asymptotic properties of the ML estima- 
tors. For small samples, they can be based on their 
approximate sampling distributions obtained through 
computer-intensive methods such as bootstrapping. 
It is usually the case that a parametric specifica- 
tion of Ao(-) is more suitable in the reliability and 
engineering situations. 

In biomedical and public health settings, it is typ- 
ical to specify Ao(-) nonparametrically, that is, to 
simply assume that Ao(-) belongs to the class of haz- 
ard rate functions with support [0, oo). This leads to 
a semiparametric model, with infinite-dimensional 
parameter Ao(-) and finite-dimensional parameters 
(a, (3, £) . Inference for this semiparametric model was 
considered in Peha, Slate and Gonzalez (2007). In 
this setting, interest is focused on the finite-dimensional 
parameters (a, (3, £) and the infinite-dimensional pa- 
rameter Ao(-) = JqXq(w) dw and its survivor func- 
tion Sq(-) = Y\' w=0 [l — Ao(dw)]. A difficulty encoun- 
tered in estimating Ao(-) is that in the intensity 
specification in (3.5), the argument of Ao(-) is the ef- 
fective age £ (s), not s, and our interest is in Ao(i), t > 
0. This poses difficulties, especially in establishing 
asymptotic properties, because the usual martin- 
gale approach as pioneered by Aalen (1978), Gill 
(1980) and Andersen and Gill (1982) (see also An- 
dersen et al., 1993 and Fleming and Harrington, 
1991) does not directly carry through. In the sim- 
ple i.i.d. renewal setting where £(s) = s — S N \r g \, 
p(k;a) = 1 and tp(w) = 1, Peha, Strawderman and 
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Hollander (2000, 2001), following ideas of Gill (1981) 
and Sellke (1988), implemented an approach using 
time transformations to obtain estimators of Ao(-) 
and So(-). In an indirect way, with partial use of 
Lenglart's inequality and Rebolledo's MCLT, they 
obtained asymptotic properties of these estimators, 
such as their consistency and their weak conver- 
gence to Gaussian processes. This approach in Peha, 
Strawderman and Hollander (2001) was also utilized 
in Peha, Slate and Gonzalez (2007) to obtain the es- 
timators of the model parameters in the more gen- 
eral model. 

The idea behind this approach is to define the 
predictable (with respect to s for fixed t) doubly in- 
dexed process t) = I{£i(s) < t},i = 1,2, ... ,n, 
which indicates whether at calendar time s the unit's 
effective age is at most t. We then define the pro- 
cesses Ni(s,t) = $°Ci{v,t)N}(dv), Ai(s,t) = 
Jq d(v,t) ■ A\(dv), and Mi(s,t) = Ni(s,t)- Ai(s,t) = 
Jq Ci(v,t)Mj (dv). Because for each t > 0, Cj(-;i) is a 
predictable and a {0, l}-valued process, then Mj(-,t) 
is a square-integrable martingale with PQV Ai{-,t). 
However, observe that for fixed s, Mj(s,-) is not 
a martingale though it still satisfies E{Mi(s,t)} = 
for every t. The next step is to have an alter- 
native expression for Ai(s,t) such that Ao(-) ap- 
pears with an argument of t instead of £i(v). With 
£ij-i(v) = £i(«)J{5y_i < v < Sij} on I{Y?(v) > 0}, 
this is achieved as follows: 

Ai(8,t;Ao(-),a,P) 

Y^v)p[Nj(v-);a} 

(4.16) • i;(Xi(v)P)I{£i(v) < t}X [£i(v)] dv 

rS, 



E 



Y?(v)p[N}(y-);a] 



j = l *'' 5 »j-l 



+ 



S t 



■ A [^-i(f)]df 
Y?(v)p[N}(v-); a ] 

■^(X i (v)P)I{S i ^ a _ ) (v)<t} 
■Ao[^t (a _)(«)]d«. 



Letting 

</Jij_i(v;a,/3) 



p[Nl{v-);a]^(X i {v)f3) 



I{Sij-i <v< Sij}, 



and defining the new "at-risk" process 
Yi(s,t;a,(3) 

= Yl I{t£(£ lj ~i(S ij - 1 +),£ lj ^ l (S ij )}} 



3=1 



(4.17) 



<Pij-i[£ijli(ty,aL,0\ 



+ I {t£(£ iN } (s _ ) (s iNHs _ ) +), 

then, after a variable transformation w = £ij-i(v) 
for each summand in (4.16), we obtain an alterna- 
tive form of Ai(s,t) given by Ai(s,t; Ao(-),a,/3) = 
/q Yi(s,w;a,P)AQ(dw). The utility of this alterna- 
tive form is that Ao(-) appears with the correct ar- 
gument for estimating it. If, for the moment, we as- 
sume that we know a and (3, by virtue of the fact 
that Mi(s, t; a, (3) has zero mean, then using the idea 
of Aalen (1978), a method-of-moments "estimator" 
of Ao(t) is 



(4.18) 



A (s*,t;a,(3) 



t I{S (s*,w;a,f3)>0} ^ ATf m 
o b (s*,w;a,p) 



8=1 



where Sq(s, t) = Y^=i Yi(s, t; a, (3). This "estimator" 
of Ao(-) can be plugged into the likelihood function 
over [0, s*] to obtain the profile likelihood of {a, (5), 
given by 

L P (s*;a,(3) 

(4.i9) = n n k? 



i=l j=l 



This profile likelihood is reminiscent of the partial 
likelihood of Cox (1972) and Andersen and Gill (1982) 
for making inference about the finite-dimensional 
parameters in the Cox proportional hazards model 
and the multiplicative intensity model. The estima- 
tor of (a, (3), denoted by (a,f3), is the maximizer of 
the mapping (a,f3) i— > Lp(s*;a, (3). Algorithms and 
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software for computing the estimate (a,P) were de- 
veloped in Peha, Slate and Gonzalez (2007). The 
estimator of Ao(i) is obtained by substituting (a,/3) 
for (a, (3) in Ao(s*, t; a, P) to yield the generalized 
Aalen-Breslow-Nelson estimator 

i / * ja [' I{S (s*,w;a,P)>0} 
Jo bo{s*,w;a,P) 

(4.20) 

n 

■J2Ni(s*,dw). 

i=l 

By invoking the product-integral representation of 
a survivor function, a generalized product-limit es- 
timator of the baseline survivor function So(t) is 
S (s*,t)=U t w = [l-Ms*,dw)]. 

Peha, Slate and Gonzalez (2007) also discussed 
the estimation of Aq(-) and the finite-dimensional 
parameters (a,P,£) in the presence of gamma-dis- 
tributed frailties. The ML estimators of these pa- 
rameters maximize the likelihood L(s*; Ao(-), a, (3, £) 
in (4.14), with the proviso that the maximizing Ao(-) 
jumps only at observed values of the £i(Sij)'s. An 
EM algorithm finds these maximizers and its imple- 
mentation is described in detail in Peha, Slate and 
Gonzalez (2007). We briefly describe the basic in- 
gredients of this algorithm here. 

With 9 = (Ao(-),a,/3,£), in the expectation step 
of the algorithm, expressions for E{Zi\6,F s *} and 
E{logZi\9,F s *}, which are easy to obtain under 
gamma frailties, are needed. For the maximization 
step, with E z \g(o) denoting expectation with respect 
to Z when the parameter vector 9 equals 9^ = 
(A[, 0) (-),a( ),/?( ),£ (0) ), we require Q{9; 0(°), F s *) = 
£ z|e( o){logL c (s*;Z,#(°)}, where L c (s;Z,9) is 

in (4.12). This Q(9;0(°\F S *) is maximized with re- 
spect to 9 = (Ao(-), a, P, £). This maximization can 
be performed in two steps: first, maximize with re- 
spect to (Ao(-),a, P) using the procedure in the case 
without frailties except that So(s,t;a, P) is replaced 
by So(s, t; Z, a, P) = Ya=x ZiYi(s, t; a, P); and second, 
maximize with respect to £ a gamma log-likelihood 
with estimated log Z{ and Z{. To start the iteration 
process, a seed value for Ao(-) is needed, which can 
be the estimate of Ao(-) with no frailties. Seed val- 
ues for (a,/3, £) are also required. Through this EM 
implementation, estimates of (Ao(-),a, are ob- 
tained and, through the product-integral represen- 
tation, an estimate of the baseline survivor function 
So(-). 



5. ILLUSTRATIVE EXAMPLES 

The applicability of these dynamic models still 
needs further and deeper investigation. We provide 
in this section illustrative examples to demonstrate 
their potential applicability. 

Example 3 . The first example deals with a data 
set in Kumar and Klefsjo (1992) which consists of 
failure times for hydraulic systems of load-haul-dump 
(LHD) machines used in mining. The data set has 
six machines with two machines each classified as 
old, medium and new. For each machine the suc- 
cessive failure times were observed and the result- 
ing data is depicted in Figure 3. Using an effec- 
tive age process £(s) = s — 5jyt( a _), this was ana- 
lyzed in Stocker (2004) (see also Stocker and Peha, 
2007) using the general class of models when the 
baseline hazard function is postulated to be a two- 
parameter Weibull hazard function Ao(i; #i, #2) = 
(i/^) 6 * 1 , and in Peha, Slate and Gonzalez (2007) 
when the baseline hazard function is nonparametri- 
cally specified. The age covariate was coded accord- 
ing to the dummy variables (Xi,X2) taking values 
(0,0) for the oldest machines, (1,0) for the medium- 
age machines and (0, 1) for the newest machines. 
The parameter estimates obtained for a nonpara- 
metric and a parametric baseline hazard function 
specification are contained in Table 1, where the 
estimates for the parametric specification are from 
Stocker (2004). The estimates of the baseline sur- 
vivor function Sq(-) under the nonpar ametric and 
parametric Weibull specifications are overlaid in Fig- 
ure 4. From this table of estimates, observe that a 
frailty component is not needed for both nonpara- 
metric and parametric specifications since the esti- 
mates of the frailty parameter £ are very large in 
both cases. Both estimates of the P\ and pz coeffi- 
cients are negative, indicating a potential improve- 
ment in the lengths of the working period of the 
machines when they are of medium age or newer, 
though an examination of the standard errors re- 
veals that we cannot conclude that the /^-coefficients 
are significantly different from zero. On the other 
hand, the estimates of a for both specifications are 
significantly greater than 1, indicating the poten- 
tial weakening effects of accumulating event occur- 
rences. From Figure 4 we also observe that the two- 
parameter Weibull appears to be a good parametric 
model for the baseline hazard function as the non- 
parametric and parametric baseline hazard function 
estimates are quite close to each other. However, a 
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Table 1 

Parameter estimates for the LHD hydraulic data for a non-parametric and a parametric 
(Weibull) specification of the baseline hazard function 



Parameter 
estimated 


With a nonparametric specification 
of A (t) 


With a parametric specification, 

A (t) = (t/0 2 ) ei 


a 


1.03 


1.03 (0.01) 


Pi 


-0.09 


-0.14 (0.20) 




-0.05 


-0.08 (0.20) 




1.54 x 10 63 


164198 (1307812) 


01 


NA 


0.97 (0.075) 


02 


NA 


0.006 (0.001) 



The values in parentheses in the third column are the approximate standard errors. 



formal procedure for validating this claim still needs 
to be developed. This is a problem in goodness-of-fit 
which is currently being pursued. 

Example 4. The second example is provided 
by fitting the general class of models to the blad- 
der cancer data in Wei, Lin and Weissfeld (1989). 
A pictorial depiction of this data set can be found 
in Peha, Slate and Gonzalez (2007), where it was 



analyzed using a nonparametric specification of the 
baseline hazard function. This data set consists of 
85 subjects and provides times to recurrence of blad- 
der cancer. The covariates included in the analysis 
are X±, the treatment indicator (1 = placebo, 2 = 
thiotepa) ; X2 , the size (in cm) of the largest initial 
tumor; and X3, the number of initial tumors. In fit- 
ting the general model in (3.5) we used p(k; a) = 



tfi - + * + + + 



5 



w - + 



4 
I 



+ + 



+ + ++ +4++ 44- 



4+ -th ft 



■+H- -H++++H- +• 



-H- ++ -H-#W- 4-Hf 44+1- + + ft 



P5 - 



+ 4-IH-4- + 4- 4- 



— h +-+H- ■+■ h-hh- + 



4--«- + + + +^ 



+ 4-* 



1000 



— i — 
2000 



3000 



400U 



calendar time 



Fig. 3. Pictorial depiction of the LHD data set from Kumar and Klefsjd (1992) which shows the successive failure occurrences 
for each of the six machines. Machines 1 and 2 have (Xi,X2) = (0,0), machines 3 and 4 have (Xi, X2) = (1,0) and machines 
5 and 6 have (Xi, X2) = (0, 1) . 
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Table 2 

Summary of estimates for the bladder data set from the Andersen-Gill (AG), Wei, Lin and Weissfeld (WLW) and Prentice, 
Williams and Peterson (PWP) methods as reported in Therneau and Grambsch (2000), together with the estimates obtained 
from the general model using two effective ages corresponding to "perfect" and "minimal" interventions 



Covariate 


Parameter 


AG 


WLW 


PWP 


General model 








marginal 


conditional 


Perfect 


Minimal 


logJV(t-) 


Q 








0.98 (0.07) 


0.79 (0.13) 


frailty 










oo 


0.97 


rx 


0i 


-0.47 (0.20) 


-0.58 (0.20) 


-0.33 (0.21) 


-0.32 (0.21) 


-0.57 (0.36) 


size 


02 


-0.04 (0.07) 


-0.05 (0.07) 


-0.01 (0.07) 


-0.02 (0.07) 


-0.03 (0.10) 


number 


03 


0.18 (0.05) 


0.21 (0.05) 


0.12 (0.05) 


0.14 (0.05) 


0.22 (0.10) 



a . Furthermore, since the data set does not con- 
tain information about the effective age, we con- 
sidered two situations for demonstration purposes: 
(i) perfect intervention is always performed [£i(s) = 
s — S iN i, a _\]j an d (h) minimal intervention is always 

performed [£i(s) = s]. The parameter estimates ob- 
tained by fitting the model with frailties, and other 
estimates using procedures discussed in the litera- 
ture, are presented in Table 2. The estimates of their 
standard errors (s.e.) are in parentheses, with the 
s.e.'s under the minimal repair intervention model 
obtained through a jackknife procedure. The other 
parameter estimates in the table are those from the 
marginal method of Wei, Lin and Weissfeld (1989), 
the conditional method of Prentice, Williams and 
Peterson (1981) and the Andersen and Gill (1982) 
method, which were mentioned earlier (cf. Therneau 



and Grambsch, 2000). Estimates of the survivor func- 
tions at the mean covariate values are in Figure 5. 

Observe from this figure that the thiotepa group 
possesses higher survivor probability estimates com- 
pared to the placebo group in either specification of 
the effective age process. Examining Table 2, note 
the important role the effective age process provides 
in reconciling differences among the estimates from 
the other methods. When the effective age process 
corresponds to perfect intervention, the resulting es- 
timates from the general model are quite close to 
those obtained from the Prentice, Williams and Pe- 
terson (1981) conditional method; whereas when a 
minimal intervention effective age is assumed, then 
the general model estimates are close to those from 
the marginal method of Wei, Lin and Weissfeld (1989). 
Thus, the differences among these existing methods 




aw eoo 



Fig. 4. Estimates of the baseline survivor function So(t) under a nonparametric and a parametric (Weibull) specification 
for the LHD hydraulic data of Kumar and Klefsjd (1992). 
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are perhaps attributable to the type of effective age 
process used. This indicates the importance of the 
effective age in modeling recurrent event data. If 
possible, it therefore behooves one to monitor and 
assess this effective process in real applications. 

6. OPEN PROBLEMS AND CONCLUDING 
REMARKS 

There are several open research issues pertaining 
to this general model for recurrent events. First is to 
ascertain asymptotic properties of the estimators of 
model parameters under the frailty model when the 
baseline hazard rate function Ao(-) is nonparamet- 
rically specified. A second problem, which arises af- 
ter fitting this general class of models, is to validate 
its appropriateness and to determine the presence of 
outlying and/or influential observations. This is cur- 
rently being performed jointly with Jonathan Quiton, 
a Ph.D. student at the University of South Carolina. 
Of particular issue is the impact of the sum-quota 
accrual scheme, leading to the issue of determin- 
ing the proper sampling distribution for assessing 
values of test statistics. This validation issue also 
leads to goodness-of-fit problems. It might, for in- 
stance, be of interest to test the hypothesis that 
the unknown baseline hazard function Ao(-) belongs 
to the Weibull class of hazard functions. In cur- 
rent research we are exploring smooth goodness-of- 
fit tests paralleling those in Peha (1998a, 1998b) and 



Agustin and Peha (2005) which build on work by 
Neyman (1937). This will lead to notions of gener- 
alized residuals from this general class of models. 
Another problem is a nonparametric Bayesian ap- 
proach to failure time modeling. Not much has been 
done for this approach in this area, though the com- 
prehensive paper of Hjort (1990) provides a solid 
contribution for the multiplicative intensity model. 
It is certainly of interest to implement this Bayesian 
paradigm for the general class of models for recur- 
rent events. 

To conclude, this article provides a selective re- 
view of recent research developments in the model- 
ing and analysis of recurrent events. A general class 
of models accounting for important facets in recur- 
rent event modeling was described. Methods of in- 
ference for this class of models were also described, 
and illustrative examples were presented. Some open 
research problems were also mentioned. 
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